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The diversity of neuron models used in contemporary theoretical neuroscience to 
investigate specific properties of covariances in the spiking activity raises the question 
how these models relate to each other. In particular it is hard to distinguish between 
generic properties of covariances and peculiarities due to the abstracted model. Here 
we present a unified view on pairwise covariances in recurrent networks in the irregular 
regime. We consider the binary neuron model, the leaky integrate-and-fire (LIF) model, 
and the Hawkes process. We show that linear approximation maps each of these models 
to either of two classes of linear rate models (LRM), including the Ornstein-Uhlenbeck 
process (OUP) as a special case. The distinction between both classes is the location of 
additive noise in the rate dynamics, which is located on the output side for spiking models 
and on the input side for the binary model. Both classes allow closed form solutions 
for the covariance. For output noise it separates into an echo term and a term due to 
correlated input. The unified framework enables us to transfer results between models. 
For example, we generalize the binary model and the Hawkes process to the situation with 
synaptic conduction delays and simplify derivations for established results. Our approach 
is applicable to general network structures and suitable for the calculation of population 
averages. The derived averages are exact for fixed out-degree network architectures and 
approximate for fixed in-degree. We demonstrate how taking into account fluctuations in 
the linearization procedure increases the accuracy of the effective theory and we explain 
the class dependent differences between covariances in the time and the frequency 
domain. Finally we show that the oscillatory instability emerging in networks of LIF models 
with delayed inhibitory feedback is a model-invariant feature: the same structure of poles 
in the complex frequency plane determines the population power spectra. 



Keywords: correlations, linear response, Hawkes process, leaky integrate-and-fire model, binary neuron, linear rate 
model, Ornstein-Uhlenbeck process 



1. INTRODUCTION 

The meaning of correlated neural activity for the processing 
and representation of information in cortical networks is still 
not understood, but evidence for a pivotal role of correlations 
increases (recently reviewed in Cohen and Kohn, 20 1 1 ) . Different 
studies have shown that correlations can either decrease (Zohary 
et al., 1994) or increase (Sompolinsky et al, 2001) the signal 
to noise ratio of population signals, depending on the readout 
mechanism. The architecture of cortical networks is dominated 
by convergent and divergent connections among the neurons 
(Braitenberg and Schtiz, 1991) causing correlated neuronal activ- 
ity by common input from shared afferent neurons in addition to 
direct connections between pairs of neurons and common exter- 
nal signals. It has been shown that correlated activity can faithfully 
propagate through convergent-divergent feed forward structures, 
such as synfire chains (Abeles, 1991; Diesmann et al., 1999), a 
potential mechanism to convey signals in the brain. Correlated 
firing was also proposed as a key to the solution of the bind- 
ing problem (von der Malsburg, 1981; Bienenstock, 1995; Singer, 
1999), an idea that has been discussed controversially (Shadlen 
and Movshon, 1999). Independent of a direct functional role 
of correlations in cortical processing, the covariance function 



between the spiking activity of a pair of neurons contains the 
information about time intervals between spikes. Changes of 
synaptic coupling, mediated by spike-timing dependent synap- 
tic plasticity (STDP, Markram et al, 1997; Bi and Poo, 1999), 
are hence sensitive to correlations. Understanding covariances in 
spiking networks is thus a prerequisite to investigate the evolution 
of synapses in plastic networks (Burkitt et al., 2007; Gilson et al., 
2009, 2010). 

On the other side, there is ubiquitous experimental evidence of 
correlated spike events in biological neural networks, going back 
to early reports on multi-unit recordings in cat auditory cortex 
(Perkel et al, 1967; Gerstein and Perkel, 1969), the observation 
of closely time-locked spikes appearing at behaviorally relevant 
points in time (Kilavik et al., 2009; Ito et al, 2011) and collec- 
tive oscillations in cortex [recently reviewed in Buzsaki and Wang 
(2012)]. 

The existing theories explaining correlated activity use a mul- 
titude of different neuron models. Hawkes (1971) developed 
the theory of covariances for linear spiking Poisson neurons 
(Hawkes processes). Ginzburg and Sompolinsky (1994) presented 
the approach of linearization to treat fluctuations around the 
point of stationary activity and to obtain the covariances for 
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FIGURE 1 | Mapping different descriptions of neuronal dynamics to 
linear rate models (LRM). The arrows indicate analytical methods which 
enable a mapping from the original spiking (LIF model, Hawkes model) or 
binary neuron dynamics to the analytically more tractable linear rate 
models. Depending on the original dynamics (spiking or binary) the 
resulting LRM contains an additive noise component x either on the output 
side (left) or on the input side (right). 



networks of non-linear binary neurons. The formal concept 
of linearization allowed Brunei and Hakim (1999) and Brunei 
(2000) to explain fast collective gamma oscillations in networks 
of spiking leaky integrate-and-fire (LIF) neurons. Correlations 
in feed-forward networks of LIF models are studied in Moreno- 
Bote and Parga (2006), exact analytical solutions for such net- 
work architectures are given in Rosenbaum and Josic (2011) 
for the case of stochastic random walk models, and thresh- 
old crossing neuron models are considered in Tchumatchenko 
et al. (2010) and Burak et al. (2009). Covariances in struc- 
tured networks are investigated for Hawkes processes (Pernice 
et al., 2011), and in linear approximation for LIF (Pernice et al., 
2012) and exponential integrate-and-fire neurons (Trousdale 
et al, 2012). The latter three works employ an expansion of the 
propagator (time evolution operator) in terms of the order of 
interaction. Finally Buice et al. (2009) investigate higher order 
cumulants of the joint activity in networks of binary model 
neurons. 

Analytical insight into a neuroscientific phenomenon based on 
correlated neuronal activity often requires a careful choice of the 
neuron model to arrive at a solvable problem. Hence a diversity of 
models has been proposed and is in use. This raises the question 
which features of covariances are generic properties of recurrent 
networks and which are specific to a certain model. Only if this 
question can be answered one can be sure that a particular result 
is not an artifact of oversimplified neuronal dynamics. Currently 
it is unclear how different neuron models relate to each other and 
whether and how results obtained with one model carry over to 
another. In this work we present a unified theoretical view on 
pairwise correlations in recurrent networks in the asynchronous 
and collective-oscillatory regime, approximating the response of 
different models to linear order. The joint treatment allows us to 
answer the question of genericness and moreover naturally leads 
to a classification of the considered models into only two cat- 
egories, as illustrated in Figure 1. The classification in addition 
enables us to extend existing theoretical results to biologically 
relevant parameters, such as synaptic delays and the presence 
of inhibition, and to derive explicit expressions for the time- 
dependent covariance functions, in quantitative agreement with 
direct simulations, which can serve as a starting point for further 
work. 

The remainder of this article is organized as follows. In the first 
part of our results in "Covariance structure of noisy rate models" 
we investigate the activity and the structure of covariance func- 
tions for two versions of linear rate models (LRM); one with input 
the other with output noise. If the activity relaxes exponentially 
after application of a short perturbation, both models coincide 
with the OUP. We mainly consider the latter case, although most 
results hold for arbitrary kernel functions. We extend the analyt- 
ical solutions for the covariances in networks of OUP (Risken, 
1996) to the neuroscientifically important case of synaptic con- 
duction delays. Solutions are derived first for general forms of 
connectivity in "Solution of the convolution equation with input 
noise" for input noise and in "Solution of convolution equa- 
tion with output noise" for output noise. After analyzing the 
spectral properties of the dynamics in the frequency domain in 
"Spectrum of the dynamics," identifying poles of the propagators 



and their relation to collective oscillations in neuronal networks, 
we show in "Population-averaged covariances'" how to obtain 
pairwise averaged covariances in homogeneous Erdos-Renyi ran- 
dom networks. We explain in detail the use of the residue theorem 
to perform the Fourier back-transformation of covariance func- 
tions to the time domain in "Fourier back transformation" for 
general connectivity and in "Explicit expression for the popula- 
tion averaged cross covariance in the time domain" for averaged 
covariance functions in random networks, which allows us to 
obtain explicit results and to discuss class dependent features of 
covariance functions. 

In the second part of our results in "Binary neurons," "Hawkes 
processes," and "Leaky integrate-and-fire neurons" we consider 
the mapping of different neuronal dynamics on either of the 
two flavors of the linear rate models discussed in the first 
part. The mapping procedure is qualitatively the same for all 
dynamics as illustrated in Figure 1: Starting from the dynamic 
equations of the respective model, we first determine the work- 
ing point described in terms of the mean activity in the net- 
work. For unstructured homogeneous random networks this 
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amounts to a mean-field description in terms of the popula- 
tion averaged activity (i.e., firing rate in spiking models). In 
the next step, a linearization of the dynamical equations is 
performed around this working point. We explain how fluc- 
tuations can be considered in the linearization procedure to 
improve its accuracy and we show how the effective linear 
dynamics maps to the LRM. We illustrate the results through- 
out by a quantitative comparison of the analytical results to 
direct numerical simulations of the original non-linear dynam- 
ics. The appendices "Implementation of noisy rate models," 
"Implementation of binary neurons in a spiking simulator 
code," and "Implementation of Hawkes neurons in a spik- 
ing simulator code." describe the model implementations and 
are modules of our long-term collaborative project to provide 
the technology for neural systems simulations (Gewaltig and 
Diesmann, 2007). 

2. COVARIANCE STRUCTURE OF NOISY RATE MODELS 
2.1. DEFINITION OF MODELS 

Let us consider a network of linear model neurons, each charac- 
terized by a continuous fluctuating rate r and connections from 
neuron j to neuron i given by the element w« of the connectivity 
matrix w. We assume that the response of neuron i to input can 
be described by a linear kernel h so that the activity in the network 
fulfills 

r(f) = h(o) * [wr(o -d) + bx(o)](t), (1) 

where f(o — d) denotes the function / shifted by the delay d, x is 
an uncorrelated noise with 

(x,(t)) = 0, {Xi(s)Xj[t)) = Syo(s - f)p 2 , (2) 

e.g., a Gaussian white noise and (f * g)(t) = f^^fit — f) 
g(t') dt' is the convolution. With the particular choice 
b = w8(o — d)* we obtain 

r(f) = [h(o) * w(r(o -d) + x(o - d))](t). (3) 

We call the dynamics (3) the linear noisy rate model (LRM) with 
noise applied to output, as the sum r + X appears on the right 
hand side. Alternatively, choosing b = 1 we define the model with 
input noise as 

r(f) = h(o) * [wr(o - d) + x(o)](f). (4) 

Hence, Equations (3) and (4) are special cases of (1). In the 
following we consider the particular case of an exponential kernel 

h(s) = -Q(s) e~ s/x , (5) 
x 

where 6 denotes the Heaviside function, 6(f) = 1 for t > 0, 0 else. 
Applying to ( 1 ) the operator O = x4- + 1 which has h as a Green's 
function (i.e., Oh = 8) we get 

d 

x — r(f) + r(f) =wr(f-d) + bx(f), (6) 
dt 



which is the equation describing a set of delay coupled Ornstein- 
Uhlenbeck-processes (OUP) with input or output noise for b = 1 
or b = w8(o — d)*, respectively. We use this representation in 
"Binary neurons" to show the correspondence to networks of 
binary neurons. 

2.2. SOLUTION OF THE CONVOLUTION EQUATION WITH INPUT NOISE 

The solution for the system with input noise obtained from the 
definition (4) after Fourier transformation is 

R = H tJ wR + HX, (7) 

where the delay is consumed by the kernel function hd(s) = 
^9(s — d)e~ ls ~ d ^ x . We use capital letters throughout the text to 
denote objects in the Fourier domain and lower case letters for 
objects in the time domain. Solved for R = (1 — Hdw) _1 HX the 
covariance function of r in the Fourier domain is found with the 
Wiener-Khinchin theorem (Gardiner, 2004) as (R(a>)R T (— w)), 
also called the cross spectrum 

C(w) = (R(w)R T (-w)> (8) 
= (1 -H d (o))w)- 1 H(w)(X( W )X T (- W )> 
H(-(o)(l - H d (-(x>)vf r r 1 

= (Hd(«) _I -wr'DCHdC-cor 1 -w T r\ 

where we introduced the matrix D = (X(w)X r (— u>)). From the 
second to the third line we used the fact that the non-delayed 
kernels H(w) can be replaced by delayed kernels Hd(d>) and that 
the corresponding phase factors e md and e~"° d cancel each other. 
If x is a vector of pairwise uncorrelated noise, D is a diagonal 
matrix and needs to be chosen accordingly in order for the cross 
spectrum (8) to coincide (neglecting non-linear effects) with the 
cross spectrum of a network of binary neurons, as described in 
"Equivalence of binary neurons and Ornstein-Uhlenbeck pro- 
cesses". 

2.3. SOLUTION OF CONVOLUTION EQUATION WITH OUTPUT NOISE 

For the system with output noise we consider the quantity y, = 
r{ + X{ as the dynamic variable representing the activity of neuron 
i and aim to determine pairwise correlations. It is easy to get from 
(3) after Fourier transformation 

R = H d w(R + X), (9) 

which can be solved for R= (1 — H l fw)~ 1 H l {wX in order to 
determine the Fourier transform of Y as 

Y = R + X=(l-H d wr 1 X. (10) 

The cross spectrum hence follows as 

C(«) = (Y( M )Y T (-oo)> (11) 
= (1 -H ( ,(oo)w)- 1 (X(co)X T (-0)))(l -HdC-c^w 7 )- 1 
= (l-HdMwr'DCl-flrff-co^r 1 , 
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FIGURE 2 | Pole structure determines dynamics. Autocovariance of the 
population activity (A,C) measured in p 2 /x and its Fourier transform called 
power spectrum (B,D) of the rate models with output noise (dots) (3) and 
input noise (diagonal crosses) (4) for delays d = 3 ms (A,B), and d = 1 ms 
(C,D). Black symbols show averages over the excitatory population activity 
and gray symbols over the inhibitory activity obtained by direct simulation. 
Light gray curves show theoretical predictions for the spectrum (20) and 
the covariance (21) for output noise and the spectrum (17) and the 
covariance (18) for input noise. Black crosses (12) in (B,D) denote the 
locations of the poles of the cross spectra - with the real parts 
corresponding to the damping (vertical axis), and the imaginary parts to 
oscillation frequencies (horizontal axis). The detailed parameters for this 
and following figures are given in "Parameters of simulations'.' 



with D = (X(eo)X r (— to)>. D is a diagonal matrix with the i-th 
diagonal entry p 2 . For the correspondence to spiking models D 
must be chosen appropriately, as discussed in "Hawkes processes" 
and "Leaky integrate-and-fire neurons" for Hawkes processes and 
LIF neurons, respectively. 

2.4. SPECTRUM OF THE DYNAMICS 

For both linear rate dynamics, with output and with input noise, 
the cross spectrum C(co) has poles at certain frequencies m 
in the complex plane. These poles are defined by the zeros of 
det(Hri(co) _1 — w) and the corresponding term with the oppo- 
site sign of 03. The zeros ofdetCH^oo) -1 -w) are solutions of the 
equation 

H^Or 1 = (1 + i<ax)e' Md = Lj 

where Lj is the j-th eigenvalue of w. The same set of poles arises 
from (1) when solving for R. For d > 0 and the exponential kernel 
(5), the poles can be expressed as 

z k (L j ) = ---W k lL j -ex\, (12) 

where Wjt is the fc-th of the infinitely many branches of the 
Lambert-W function (Corless et al., 1996). For vanishing synaptic 
delay d = 0 there is obviously only one solution for every Lj given 
byz=f(L i -l). 

Given the same parameters d, w, t, the pole structures of 
the cross spectra of both systems (8) and (11) are identical, 
since the former can be obtained from the latter by multiplica- 
tion with (Hd(a))Hrf(-co)) _1 = (H(co)H(-co)) -1 , which has no 
poles. The only exception causing a different pole structure for the 
two models is the existence of an eigenvalue Lj = 0 of the connec- 
tivity matrix w, corresponding to a pole z(0) = i . However, this 
pole corresponds to an exponential decay of the covariance for 
input noise in the time domain and hence does not contribute to 
oscillations. For output noise, the multiplication with the term 
(H(w)H(— w)) -1 , vanishing at oj = K cancels this pole in the 
covariance. Consequently both dynamics exhibit similar oscilla- 
tions. A typical spectrum of poles for a negative eigenvalue Lj < 0 
is shown in Figures 2B,D. 

2.5. POPULATION-AVERAGED C0VARIANCES 

Often it is desirable to consider not the whole covariance matrix 
but averages over subpopulations of pairs of neurons. For instance 
the average over the whole network would result in a single scalar 
value. Separately averaging pairs, distinguishing excitatory and 
inhibitory neuron populations, yields a 2 by 2 matrix of covari- 
ances. For these simpler objects closed form solutions can be 
obtained, which already preserve some useful information and 
show important features of the network. Averaged covariances 
are also useful for comparison with simulations and experimental 
results. 

In the following we consider a recurrent random network of 
N e = N excitatory and N{ = yN inhibitory neurons with synap- 
tic weight w for excitatory and — gw for inhibitory synapses. The 
probability p determines the existence of a connection between 



two randomly chosen neurons. We study the dynamics aver- 
aged over the two subpopulations by introducing the quantities 

r « = W a E;efl r ; andn0isetermSX a = W a T,j€a X j{° Ia € 

indices 2 and £ stand for inhibitory and excitatory neurons 
and corresponding quantities. Calculating the average local input 
IZjea w jk r k to a neuron of type a, we obtain 

Na 1 E E w m = N a' E E w m + E E w m (13) 

j€a k \j€ak<E£ jea keX ) 

= N a X ipNaW E n ~ pN a gW E r k I 
\ ke£ keX I 

= pwN(r £ - ygrz), 

where, from the second to the third line we used the fact that 
in expectation a given neuron k has pN a targets in the popula- 
tion a. The reduction to the averaged system in (13) is exact if 
in every column k in w» there are exactly K non-zero elements 
for j € £ and yK for j e X, which is the case for networks with 
fixed out-degree (number of outgoing connections of a neuron to 
the neurons of a particular type is kept constant), as noted earlier 
(Tetzlaff et al., 2012). For fixed in-degree (number of connec- 
tions to a neuron coming in from the neurons of a particular type 
is kept constant) the substitution of rj €a by r a is an additional 
approximation, which could be considered as an average over pos- 
sible realizations of the random connectivity. In both cases the 
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FIGURE 3 | Limits of the theory for fixed in-degree and fixed 
out-degree. Autocovariance (A) and covariance (B) in random networks 
with fixed in-degree (dots) and fixed out-degree (crosses). Simulation 
results for cee, C£x. and Cn are shown in dark gray, black and light gray, 
respectively for synaptic weight w = 0.011 far from bifurcation. For larger 
synaptic weight w = 0.018 close to bifurcation (see text at the end of 
"Population-averaged covariances"), cee is also shown in (D) for fixed 
in-degree (dark gray dots) and for fixed out-degree (black dots). 
Corresponding theoretical predictions for the autocovariance (34) (A) and 
the covariance (18) (B,D) are plotted as light gray curves throughout. The 
set of eigenvalues is shown as black dots in panel (C) for the smaller 
weight. The gray circle denotes the spectral radius vv v /Wp(1 — p)(1 + yg 2 ) 
(Rajan and Abbott, 2006; Kriener et al., 2008) confining the set of 
eigenvalues for the larger weight. The small filled gray circle and the 
triangle show the effective eigenvalues L of the averaged systems for small 
and large weight, respectively. 



effective population-averaged connectivity matrix M turns out 
to be 

M = Kw( 1 ~ VS \ (14) 

V 1 -yg) 

with K = pN. So the averaged activities fulfill the same Equations 
(3) and (4) with the non-averaged quantities r, x, and w replaced 
by their averaged counterparts f = (rg, rj) T , x = (xe, xj) T , and 
M. The population averaged activities r a are directly related to 

) , with 

CIS CTT ) 

c ab = N- l N l ; 1 j: i€a j: }eb c lj .mth 

Dab = N- 1 N; l Ij2x l J2^\ (15) 

\ i € a j€b I 

i e a j e b 
= KbN a /N 2 a p 2 = h ab N- 1 p 2 

we replace D by D = p 2 ( ^ " , | and c by c so that the 

same Equations (11) and (8) and their general solutions also hold 
for the block-wise averaged covariance matrices. 

The covariance matrices separately averaged over pairs of 
excitatory, inhibitory or mixed pairs are shown in Figure 2 for 
both linear rate dynamics (3) and (4). (Parameters for all sim- 
ulations presented in this article are collected in "Parameters 
of simulations," the implementation of LRM is described in 
"Implementation of noisy rate models"). The poles of both mod- 
els shown in Figure 2B are given by (12) and coincide with the 
peaks in the cross spectra (8) and ( 1 1 ) for output and input noise, 
respectively. The results of direct simulation and the theoretical 
prediction are shown for two different delays, with the longer 
delay leading to stronger oscillations. 

Figure 3C shows the distribution of eigenvalues in the com- 
plex plane for two random connectivity matrices with different 
synaptic amplitudes w. The model exhibits a bifurcation, if at least 
one eigenvalue assumes a zero real part. For fixed out-degree the 
averaging procedure (13) is exact, reflected by the precise agree- 
ment of theory and simulation in Figure 3D. For fixed in-degree, 
the averaging procedure (13) is an approximation, which is good 
only for parameters far from the bifurcation. Even in this regime 
still small deviations of the theory from the simulation results are 
visible in Figure 3B. On the stable side close to a bifurcation, the 
appearance of long living modes causes large fluctuations. These 
weakly damped modes appearing in one particular realization of 
the connectivity matrix are not represented after the replacement 
of the full matrix w by the average M over matrix realizations. The 
eigenvalue spectrum of the connectivity matrix provides an alter- 
native way to understand the deviations. By the averaging the set 
of N eigenvalues of the connectivity matrix is replaced withby the 
two eigenvalues of the reduced matrix M, one of which is zero due 
to identical rows of M. The eigenvalue spectrum of the full matrix 



is illustrated in Figure 3C. Even if the eigenvalue(s) L of M are 
far in the stable region (corresponding to 3(z(L M )) > 0) some 
eigenvalues L w of the full connectivity matrix in the vicinity of 
the bifurcation region may still have an imaginary part becom- 
ing negative and the system can feel their influence, shown in 
Figure 3D. 

2.6. FOURIER BACK TRANSFORMATION 

Although the cross spectral matrices (8) and (11) for both dynam- 
ics look similar in the Fourier domain, the procedures for back 
transformation differ in detail. In both cases, the Fourier inte- 
gral along the real w-axis can be extended to a closed integra- 
tion contour by a semi-circle with infinite radius centered at 
0 in the appropriately chosen half-plane. The half-plane needs 
to be selected such that the contribution of the integration 
along the semi-circle vanishes. By employing the residue theorem 
(Bronstein et al, 1999) the integral can be replaced by a sum over 
residua of the poles encircled by the contour. For a general covari- 
ance matrix we only need to calculate c(f) for f > 0, as for t < 0 
the solution can be found by symmetry c(f) = c r (— f). 

For input noise it is possible to close the contour in the upper 
half-plane where the integrand C(m) e' wt vanishes for |co| — > oo 
forallf > 0,as |C«(ci))| decays as |mp 2 . This can be seen from (8), 
because the highest order of HJ 1 oc a> appearing in de^ffj 1 — 
w) is equal to the dimensionality N of w (N = 2 for M), and in 
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det(adjugate matrix ij of HJ — w) it is N — 1 (i = j) or N — 2 
So KHJ 1 -w) _1 | is proportional to |to| _1 |e _ia>d | and 

|C(co)| oc |cop 2 for large | co| . 

For the case of output noise (11) C(co) can be obtained 
from the C(co) for input noise (8) multiplied with 
(Hd(a>)H d (— co)) -1 ~ | co | 2 for large |co|. The multiplication 
with this factor changes the asymptotic behavior of the inte- 
grand, which therefore contains terms converging to a constant 
value and terms decaying like Icop 1 for |co| — > oo. These terms 
result in non-vanishing integrals over the semicircle in the upper 
half-plane and have to be considered separately. To this end we 
rewrite (11) as 



C(o>) = ((1 - H d (a>)w)- 1 H d ((x>)w+ 1) 

D(w r H d (-co)(l - HdC-c^w 7 )- 1 + 1) (16) 
= (1 -H rf (a ) )w)- 1 H rf (o))wDw r H lJ (-oo)(l -H d (-oo)w r )- 1 
+ (1 - H d ( M )w)- 1 H d (oo)wD 
+ Dw r H rf (-co)(l - HrfC-^w 7 )" 1 
+ D, 



and find the constant term D which turns into a 8-function in the 
time domain. The first term in the second line of (16) decays like 
I cop 2 and can be transformed just as C(co) for input noise closing 
the contour in the upper half-plane. The second and third term 
are the transposed complex conjugates of each other, because of 
the dependence of H on — co instead of co, and require a special 
consideration. Multiplied by e"*' under the Fourier integral, the 
first term is proportional to H^e""' ~ co _1 e* w< - t_ '^ and vanishes 
faster than | to | 1 for large |oj| in the upper half-plane for t > d 
and in the lower half plane for t < d. For the second term the half 
planes are interchanged. The application of the residue theorem 
requires closing the integration contour in the half-plane where 
the integral over the semi-circle vanishes faster than loop 1 . For 
w = M and in the general case of a stable dynamics all poles of 
the first term are in the upper half-plane Zs(Zk(Lj)) > 0, and have 
no contribution to c(f) for t < d. For the second term the same is 
true for t > —d; these terms correspond to the jumps of c(t) after 
one delay, caused by the effect of the sending neuron arriving at 
the other neurons in the system after one synaptic delay. These 
terms correspond to the response of the system to the impulse 
of the sending neuron - hence we call them "echo terms" in the 
following (Helias et al, 2013). The presence of such discontinu- 
ous jumps at time points d and — d in the case of output noise is 
reflected in the convolution of kw with D in the time domain in 
(37). For input noise the absence of discontinuities can be inferred 
from the absence of such terms in (33), where the derivative of 
the correlation function is equal to the sum of finite terms. The 
first summand in (16) corresponds to the covariance evoked by 
fluctuations propagating through the system originating from the 
same neuron and we call it "correlated input term". In the system 
with input noise a similar separation into effective echo and cor- 
related input terms can be performed. We obtain the correlated 
input term as the covariance in an auxiliary population without 
outgoing connections and echo terms as the difference between 



the full covariance between neurons within the network and the 
correlated input term. 

2.7. EXPLICIT EXPRESSION FOR THE POPULATION AVERAGED CROSS 
COVARIANCE IN THE TIME DOMAIN 

We obtain the population averaged cross spectrum in a recurrent 
random network of Ornstein-Uhlenbeck processes (OUP) with 
input noise by inserting the averaged connectivity matrix w = M 
(14) into (8). The explicit expression for the covariance func- 
tion follows by taking into account all (both) eigenvalues of M 
with values 0 and L = Kw(l — yg). The detailed derivation of the 
results presented in this section are documented in "Calculation 
of the Population Averaged Cross Covariance in Time Domain". 
The expression for the cross spectrum (8) takes the form 



CO) 



(17) 



where we introduced /(co) = (H^cop 1 — Lp 1 as a short 
hand. Sorting the terms by their dependence on co, intro- 
ducing the functions <t>i(co), . . . , <t>4(co) for this dependence, 
and cpi(t), . . . , cp,j(f) for the corresponding functions in the 
time domain, the covariance in the time domain c(t) = 



1 r+oo 



C(co)e ,mf dco takes the form 



c(t) = Dcpi(f) 



The previous expression is valid for arbitrary D. In simulations 
presented in this article we consider identical marginal input 
statistics for all neurons. In this case the averaged activities for 
excitatory and inhibitory neurons are the same, so we can insert 
the special form of D given in (15), which results in 



c(f) 



£ ( 1 0 

N \0 y- 1 



<Pi(0 



-I Kw 

N 



[yg -g \ 



(18) 



N \-g-Y 



N \ g y 



The time-dependent functions cpi,...,cp4 are the same 
in both cases. Using the residue theorem cp,(f) = 
^/-^ $ .( w ) e '' Ut ^ = i Ez ep olesof*,Res(<I',,z) e ' z ' for f^O 
they can be expressed as a sum over the poles zjt(L) given by (12) 
and the pole z = ^ of H^(co). At co = z^(I) the residue of/(co) is 

Res(f, co = Zk(L)) = (idL + ixe"" d y l , the residue of Hrf(oo) at 
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FIGURE 4 | Functions cp 0 (D), <p-, (C), cp 2 , <p 3 (B), <p 4 (A) introduced in 
(19) and (22) for decomposition of covariance c(f). In (B) <P3 (— 0 is 

shown in gray and <P2(0 in black. The two functions are continuations of 
each other, joint at f = 0. Both functions appear in the echo term for 
input noise. The function <po in (D) describing the corresponding echo 
term in the case of output noise is shifted to be aligned with the 
function in (B) to facilitate the comparison of (B,D). Parameters in all 
panels are d = 3 ms, x = 10 ms, L = —1 .72. 



z = | is — ^e d / x , so that the explicit forms of (pi 94 follow as 

(pi(0 = !'Res(/, M )/(-oj)e' w ' 

m=Zfc(£) 

<P2(t) = ! ' Res (/' <*>)/(- w)H d ( M )e ,W 

a>=Zk(L) 

+^'(0'(-0 

93(f) = J] iRes(f, co)/(-o)if d (-co)e io,t (19) 

<P4(t) = ! ' Res (/' M)/(-w)H d (to)H d (- W )e ,Mt 

o)=z*(i) 

♦^(0'(-0- 

The corresponding expression for C(ou) for output noise is 
obtained by multiplying (17) with HJ 1 (utyHJ 1 (—a>) = (1 + 

W 2 T 2 ) 

C(u>) = HJ^)HJ \-co)/(oo)/(-co) (20) 
x (1 + Kw ( Y ^ ~Y ) H d(«>))D(l + «v f ™ g _\ J H rf (-»)), 

which, after Fourier transform, provides the expression for c(f) in 
the time domain for t ^ 0 

c(f) = MDM T (pi(f) + MD<p 0 (f) + D8(f) 

4 (:;.>«>■ 

As in (18), the first line holds for arbitrary D, and the second for 
D given by (15), valid if the firing rates are homogeneous, 91 is 
defined as before, and 

y 0 (t) = Q(t-d) (dL + xe^y 1 e ,Mt (22) 

G) = Zjt(£) 

vanishes for t < d. All matrix elements of the first term in (21) are 
identical. Therefore all elements of c(f) are equal for 0 < \t\ < d. 
Both rows of the matrix in front of cpo are identical, so for t > 0 
the off diagonal term cxs coincides with egg and cgx with cxx 
and vice versa for t < 0. 

As an illustration we show the functions cpo, . . . , 94 for one 
set of parameters in Figure 4. The left panels (A,C) correspond 
to contributions to the covariance caused by common input to a 
pair of neurons, the right panels (B,D) to terms due to the effect 
of one of the neurons' activities on the remaining network (echo 
terms). The upper panels (A,B) belong to the model with input 
noise, the lower panel (C,D) to the one with output noise. 



For the rate dynamics with output noise, the term with 91 in 
(21) (shown in Figure 4C) is symmetric and describes the com- 
mon input covariance and the term with 90 (shown in Figure 4D) 
is the echo part of the covariance. For the rate dynamics with 
input noise (18) the term containing 94 (shown in Figure 4A) is 
caused by common input and is hence also symmetric, the terms 
with 92 and 93 (shown in Figure 4B) correspond to the echo 
part and have hence their peak outside the origin. The sec- 
ond echo term in (18) is equal to the first one transposed and 
with opposite sign of the time argument, so we show 92(f) 
and 93 (— f) together in one panel in Figure 4B. Note that for 
input noise, the term with 91 describes the autocovariance, which 
corresponds to the term with the o-function in case of output 
noise. 

The solution (18) is visualized in Figure 6, the solution (21) 
in Figure 7, and the decomposition into common input and 
echo parts is also shown and compared to direct simulations in 
Figure 8. 

3. BINARY NEURONS 

In the following sections we study, in turn, the binary neuron 
model, the Hawkes model and the LIF model and show how they 
can be mapped to one of the two OUPs; either the one with input 
or the one with output noise, so that the explicit solutions (18) 
and (21) for the covariances derived in the previous section can 
be applied. In the present section, we start with the binary neuron 
model (Ginzburg and Sompolinsky, 1994; Buice et al., 2009). 

Following Ginzburg and Sompolinsky (1994) the state of the 
network of N binary model neurons is described by a binary 
vector n e {0, 1} N and each neuron is updated at independently 
drawn time points with exponentially distributed intervals of 
mean duration x. This stochastic update constitutes a source 
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of noise in the system. Given the z'-th neuron is updated, the 
probability to end in the up-state (n; = 1) is determined by 
the gain function -F,(n) which depends on the activity n of 
all other neurons. The probability to end in the down state 
(m = 0) is 1 — -F,( n )- Here we implemented the binary model 
in the NEST simulator (Gewaltig and Diesmann, 2007) as 
described in "Implementation of Binary Neurons in a Spiking 
Simulator Code". Such systems have been considered earlier 
(Ginzburg and Sompolinsky, 1994; Buice et al., 2009), and here 
we follow the notation employed in the latter work. In the fol- 
lowing we collect results that have been derived in these works 
and refer the reader to these publications for the details of the 
derivations. The zero-time lag covariance function is defined 
as qj(t) = (ni(t)tij(t)) — ai(t)aj(t), with the expectation value () 
taken over different realizations of the stochastic dynamics. Here 
a(f) = («i(f), • • ■ , fljv(f)) T is the vector of mean activities a,(f) = 
(«;(t)). Cij(t) fulfills the differential equation 



X Jt C " (f) = ~ 2c, i (t) + {(n ) (t) ~ ^W^W) 
+ ((m(t)-ai(t))Fj(n)). 

In the stationary state, the covariance therefore fulfills 

Cij = - ((«; - aj)F,(n)} + - {(«; - ai)Fj(n)). (23) 

The time lagged covariance Cy (f, s) = («,(t)»;(s)) — aj(t)dj(s) ful- 
fills for t > s the differential equation 



•t-Cy(f, s) = -dj{t, s) + (Fi(n, tKnj(s) - a ; (s))>. (24) 

This equation is also true for i = j, the autocovariance. The term 
(F,(n, t)(nj(s) — flj(s))) has a simple interpretation: it measures 
the influence of a fluctuation of neuron j at time s around its 
mean value on the gain of neuron z at time t (Ginzburg and 
Sompolinsky, 1994). We now assume a particular form for the 
coupling between neurons 



Fi(n, t) = cKJ«n(f - d)) 



N 
k= 1 



(25) 



where J, is the vector of incoming synaptic weights into neuron 
i and <p is a non-linear gain function. Assuming that the fluctu- 
ations of the total input J,n into the z'-th neuron are sufficiently 
small to allow a linearization of the gain function (f>, we obtain 
the Taylor expansion 

Fi(n, t) = F,(a) + c)/(J,a) J«(n(f - d) - a(f - d)), 

where 

*'(J,a) (26) 
is the slope of the gain function at the point of mean input. 



Up to this point the treatment of the system is identical to 
the work of Ginzburg and Sompolinsky (1994). Now we present 
an alternative approach for the linearization which takes into 
account the effect of fluctuations in the input. For sufficiently 
asynchronous network states, the fluctuations in the input J,n(f — 
d) to neuron z can be approximated by a Gaussian distribution 
Af(\i, a). In the following we consider a homogeneous ran- 
dom network with fixed in-degree as described in "Population- 
averaged covariances". As each neuron receives the same number 
K of excitatory and yK inhibitory synapses, the marginal statis- 
tics of the summed input to each neuron is identical. The mean 
input to a neuron then is \i = KJ(l — yg)a, where a is the mean 
activity of a neuron in the network. If correlations are small, the 
variance of this input signal distribution can be approximated as 
the sum of the variances of the individual contributions from the 
incoming signals, resulting in a 2 = KJ 2 (l + yg 2 ) a{\ — a), where 
we used the fact that the variance of a binary variable with mean 
a is a(l — a). This results from a direct calculation: since n S 
{0, 1}, n 2 = n, so that the variance is (n 2 ) — (n) 2 = (n) — (n) 2 = 
a(l — a). Averaging the slope <\>' of the gain function over the 
distribution of the input variable results in the averaged slope 



/oo 
JV([L,0, X) ((/(*) 
-oo 

with a, x) ■■ 



dx 



(27) 



1 



27TCT 



exp 



2a 2 



The two alternative methods of linearization of c(> are illustrated in 
Figure 5. In the given example, the linearization procedure tak- 
ing into account the fluctuations of the input signal results in a 
smaller effective slope (<\>'} than taking the slope ())'(«) at the mean 
activity a near its maximum. Averaging the slope (rj/) over this 
distribution fits simulation results better than 4>'(a) calculated at 
the mean of a, as shown in Figure 6. 

The finite slope of the non-linear gain function can be under- 
stood as resulting from the combination of a hard threshold with 
an intrinsic local source of noise. The inverse strength of this noise 




FIGURE 5 | Alternative linearizations of the binary neuron model. The 

black curve represents the non-linear gain function <|>(x) = j + j tanh(fix). 
The dashed gray line is its tangent at the mean input value (denoted by the 
diagonal cross). The solid curve is the slope (<|>') averaged over the 
distribution of the fluctuating input (27). This distribution estimated from 
direct simulation is presented by black dots, the corresponding theoretical 
prediction of a normal distribution Af(\i, a) (27) is shown as the light gray 
curve. 
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FIGURE 6 | Binary model neuron corresponds to OUP model with 
input noise. Autocovariance (A), crosscovarince (B), and autocovariance of 
population averaged activity (C,D), for binary neurons (dots) and rate model 
with input noise (crosses), cee, C£x and cxx are shown in black, gray, and 
light gray. Corresponding theoretical predictions (18) in (CD), (34) in (A), 
their difference in (C) are plotted as light gray curves throughout. Dashed 
curve in (C) represents the theoretical prediction using the linearization 
with the slope at the mean activity (26), the solid curve shows the results 
for the slope averaged over Gaussian distributed input fluctuations (27). The 
spread of the simulation results for binary neurons in (C) is due to different 
realizations of the random connectivity. (E,F) are the same as (A,B) but for 
the presence of a synaptic delay d = 10 ms instead of d = 0.1 ms. 



determines the slope parameter P (Ginzburg and Sompolinsky, 
1994). In this sense, the network model contains two sources of 
noise, the explicit local noise, quantified by P and the fluctuating 
synaptic input interpreted as self-generated noise on the network 
level, quantified by cr. Even in the absence of local noise (P — > oo), 
the above mentioned linearization is applicable and yields a finite 
effective slope {§') (27). In the latter case the resulting effective 
synaptic weight is independent of the original synapse strength 
(Grytskyy etal, 2013). 

We now extend the classical treatment of covariances in binary 
networks (Ginzburg and Sompolinsky, 1994) by synaptic con- 
duction delays. In (25) F,(n, f) must therefore be understood as 
a functional acting on the function n(f') for t' e [— oo, f], so 
that also synaptic connections with time delay d can be realized. 
We define an effective weight vector to absorb the gain factor as 
w,- = P,J,, with either p, = cp'(|x) or P; = depending on the 
linearization procedure, and expand the right hand side of (24) to 
obtain 

N 

(f;(n, t)(nj(s) - cij(s))) = ^2 WikCkjit - d, s). 

k= 1 



Thus the cross-covariance fulfills the matrix delay differential 
equation 

d 

x— c(t, s) + c(f, s) = wc(f - d, s). (28) 
dt 

This differential equation is valid for t > s. For the stationary 
solution, the differential equation only depends on the relative 
timing u = t — s 

d 

x— c(w) + c(m) =wc(u-d). (29) 
du 

The same linearization applied to (23) results in the boundary 
condition for the solution of the previous equation 

2c(0) = wc(-d) + (wc(-d)) T (30) 

or, if we split c into its diagonal and its off-diagonal parts c a 
and c^l 

2c^(0) = wc^(-d) + (wc^{-d)) T + O (31) 
with O = wc a (-d) + (wc a (-d)) T . 

In the following section we use this representation to demonstrate 
the equivalence of the covariance structure of binary networks to 
the solution for OUP with input noise. 

3.1. EQUIVALENCE OF BINARY NEURONS AND 
ORNSTEIN-UHLENBECK PROCESSES 

In the following subsection we show that the same Equations 
(29) and (31) for binary neurons also hold for the Ornstein- 
Uhlenbeck process (OUP) with input noise. In doing so here 
we also extend the existing framework of OUP (Risken, 1996) 
to synaptic conduction delays d. A network of such processes is 
described by 

d 

x—r(t) + r(t)=wt(t-d)+x(t), (32) 
at 

where x is a vector of pairwise uncorrelated white noise with 
(x(f))x = 0 and (x,(t)xj(t + t')) x = o y 8(f')p 2 . With the help of 
the Green's function G satisfying (t,^ + 1) G(f) = 5(f), namely 
G(r) = i 9(f) e~ f / T , we obtain the solution of Equation (32) as 

r(f) = tG(f)r(0) + f G(t- f')(wr(f' - d) + x(f')) dt' . 
Jo 

The equation for the fluctuations 8r(f) = r(f) — (r(f)) x around 
the expectation value 

8r(f) = f G(f - f')(w8r(f' - d) + x(f')) dt' 
Jo 

coincides with the noisy rate model with input noise (4) with 
delay d and convolution kernel h = G. In the next step we inves- 
tigate the covariance matrix c«(f, s) = (8r,(f + s)hrj(f)) x to show 
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for which choice of parameters the covariance matrices for the 
binary model and the OUP with input noise coincide. To this end 
we derive the differential equation with respect to the time lag s 
for positive lags s > 0 

thefts) = |t^5r(f + s)8r r (f)| x (33) 

= ((w8r(f + s-d)- Sr(f + 5) + x(f + s))Sr r (f))x 
= wc(t, s — d) — c(t, s), 

where we used (x(f + s))8r(f)) x = 0, because the noise is real- 
ized independently for each time step and the system is causal. 
Equation (33) is identical to the differential equation satisfied 
by the covariance matrix (28) for binary neurons (Ginzburg and 
Sompolinsky, 1994). To determine the initial condition of (33) we 
need to take the limit c(f , 0) = lim s _^ +u c(f , s). This initial condi- 
tion can be obtained as the stationary solution of the following 
differential equation 

T— c(f,0)= Km (It — 8r(f + s)8r r (f) ) + ( Sr(t + s)x — 5r r (f) ) | 
dt ^+o\\ dt l x \ dt I J 

= Km ( ((wSrft + s-d)- Sr(t + s) + x(t + s))Sr r (f))x 

s->+0 V 

+(8r(f + s)(8r r (t - d)w T - hr T (t) + x r (t))>x) 
= -2c(f, 0) + wc(f, -d) + c(f - d, d)w T + D. 

Here we used that (x(f + s)8r r (f)) vanishes due to independent 
noise realizations and causality and 

D = Km (8r(t + s)x T (t)) x 

s->+0 

rt+s 

= Km / G(f + s- f')(w(8r(f' - d)x T (t)) x + (x{t')x T (t)) x )dt' 

5^+0, s<dj 0 • , < 1 , ' 

=0 causality =15(f-t')p 2 

rt+s 

= Km / G(t + s- t')U(t- t')p 2 dt' 

s<dj 0 

= Km G(s)lp 2 = -lp 2 . 

s<d T 

In the stationary state, c only depends on the time lag s and is 
independent of the first time argument t , which, with the symme- 
try c(— d) T = c(d) yields the additional condition for the solution 

of (33) 

2c(0) = wc(-cZ) + (wc(-d)) T + D 

or, if c is split in diagonal and off-diagonal parts c a and c^, 
respectively, 

2c^(0) = wc^(-d) + (wc^{-d)) T + O 
2c fl (0) = wc^{-d) + (wc^(-d)f + D 

with O = wc„(— d) + (wc fl (— d)) T . In the equation for the auto- 
covariance c fl the first two terms are contributions due to the cross 
covariance. In the state of asynchronous network activity with 



Cjj ~ for i ytz j these terms are typically negligible in compar- 
ison to the third term because ^ k WfcCjd ~ \vKN~ 1 = pw, which 
is typically smaller than 1 for small effective weights w < 1 and 
small connection probabilities p 1. In this approximation with 
(33) the temporal shape of the autocovariance function is expo- 
nentially decaying with time constant t. With c a (0) ~ D/2 the 
approximate solution for the autocovariance is 

c a (t) = -expf-^J. (34) 

The cross covariance then satisfies the initial condition 

2c^(0) = wc^i-d) + (wc^i-d)) 7 + O 
O = wD/2+ (wD/2) T , 

which coincides with (31) for binary neurons if the diagonal 
matrix containing the zero time autocorrelations c fl (0) for binary 
neurons is equal to D/2, i.e., if the amplitude of the input 
noise p 2 = 2t«(l — a) and the effective linear coupling satisfies 
w; = PJ,-. Figure 6 shows simulation results for population aver- 
aged covariance functions in binary networks and in networks of 
OUPs with input noise where the parameters of the OUP net- 
work are chosen according to the requirements derived above. 
The theoretical results (18) agree well with the direct simulations 
of both systems. For comparison, both methods of linearization, 
as explained above, are shown. The linearization procedure which 
takes into account the noise on the input side of the non-linear 
gain function results in a more accurate prediction. Moreover, 
the results derived here extend the classical theory (Ginzburg and 
Sompolinsky, 1994) by considering synaptic conduction delays. 
Figure 8 shows the decomposition of the covariance structure for 
a non-zero delay d = 3 ms. For details of the implementation see 
"Implementation of binary neurons in a spiking simulator code". 
The explicit effect of introducing delays into the system, such as 
the appearance of oscillations in the time dependent covariance, is 
presented in (E,F) of Figure 6, differing from (A,B) of this figure, 
respectively, only in the delay (d = 10 ms for (E,F), d = 0.1 ms 
for (A,B))- 

4. HAWKES PROCESSES 

In the following section we show that to linear order the covari- 
ance functions in networks of Hawkes processes (Hawkes, 1971) 
are equivalent to those in the linear rate network with output 
noise. Hawkes processes generate spikes randomly with a time 
density given by r(f), where neuron i generates spikes at a rate 
r,(t), realized independently within each infinitesimal time step. 
Arriving spike trains s influence r according to 

r(f) = V+(h d *Js)(t), (35) 

with the connectivity matrix J and the kernel function hd includ- 
ing the delay. Here v is a constant base rate of spike emis- 
sion assumed to be equal for each neuron. Here we employ 
the implementation of the Hawkes model in the NEST simu- 
lator (Gewaltig and Diesmann, 2007). The implementation is 
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described in "Implementation of Hawkes neurons in a spiking 
simulator code". 

Given neuron j spiked at time u < t, the probability of a spike 
in the interval [f, t + 8f) for neuron i is 1 if i = j, u = t (the 
neuron spikes synchronously with itself) and r,(f)8f + o(8f 2 ) oth- 
erwise. Considering the system in the stationary state with the 
time averaged activity r = (s(f)) we obtain a convolution equa- 
tion for time lags x > 0 for the covariance matrix with the entry 
Cij(x) for the covariance between spike trains of neurons i and; 

c(t) = (s(f + T)s r (f)> - (s(f + t)>(s T (f)> (36) 
= ((o(t)l + r(f-|-T:))s T (f)> -?f T 
= (r(t + t)(s r (f)-? T )> + D F 
= {(y+(h d * Js)(f + t))(s r (f) - f T )) + D f 

= h d *](s(t+x)(s T (t)-i T ))+Dr 

= (h d *Jc)(x) + D f , 

with the diagonal matrix D F = 8(t)diag(r), which has been 
derived earlier (Hawkes, 1971). If the rates of all neurons are 
equal, r, = r, all entries in the diagonal matrix are the same, 
Df = 8(t)lr. In the subsequent section we demonstrate that the 
same convolution Equation (36) holds for the linear rate with 
output noise. 

4.1. CONVOLUTION EQUATION FOR LINEAR NOISY RATE NEURONS 

For the linear rate model with output noise we use Equation 
(3) for time lags x > 0 to obtain a convolution equation for the 
covariance matrix of the output signal vector y = r + x as 



c(t) = (y(t + t)(y T (f)-f T )> 



(37) 



= {Qi d * wy + x)(t + t)(y r (f) - r r )> 

= Qi d * wc)(t) + (x(f + t)(r T (f) - f T )> + (x(f + t)x T (f)} 

= (ft d *wc)(t) + D, 

where we utilized that due to causality the random noise sig- 
nal generated at t + x has no influence on r(f), so the respective 
correlation vanishes. D is the covariance of the noise as in (11), 
Dijix) = (x t (t)Xj(t + t)) = 8y8(t)p 2 . If p is chosen such that p 2 
coincides with the averaged activity r in a network of Hawkes neu- 
rons and the connection matrix w is identical to J of the Hawkes 
network, the Equations (36) and (37) are identical. Therefore the 
cross spectrum of both systems is given by (11). 

4.2. NON-LINEAR SELF-CONSISTENT RATE IN RECTIFYING HAWKES 
NETWORKS 

The convolution Equation (36) for the covariance matrix of 
Hawkes neurons is exact if no element of r is negative, which 
is particularly the case for a network of only excitatory neu- 
rons. Especially in networks including inhibitory couplings, the 
intensity r; of neuron i may assume negative values. A neu- 
ron with r, < 0 does not emit spikes, so the instantaneous 
rate is given by X; = [r,(f)]+ = 6(n(0) r,-(f), with the Heaviside 
function 9. We now take into account this effective nonlinearity 



-the rectification of the Hawkes model neuron- in a similar man- 
ner as we already used to linearize binary neurons. If the network 
is in the regime of low spike rates, the fluctuations in the input 
of each neuron due to the Poissonian arrival of spikes are large 
compared to the fluctuations due to the time varying intensities 
r(f). Considering the same homogeneous network structure as 
described in "Population-averaged covariances," the input statis- 
tics is identical for each cell i, so the mean activity Xo = (X,-) 
is the same for all neurons i. The superposition of the synap- 
tic inputs to neuron i cause an instantaneous intensity r; that 
follows approximately a Gaussian distribution J\f([i, a, r;) with 
mean \i = (r) = v + Xo_fC/(l — gy) and standard deviation a = 



1 A/ It ^(1 + 8 2 y)- These expressions hold for the 



n - m 

exponential kernel (5) due to Campbell's theorem (Papoulis and 
Pillai, 2002), because of the stochastic Poisson-like arrival of 
incoming spikes, where the standard deviation of the spike count 
is proportional to the square root of the intensity X.rj. The rate Xo 
is accessible by explicit integration over the Gaussian probability 
density as 
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This equation needs to be solved self-consistently (numerically 
or graphically) to determine the rate in the network, as the 
right hand side depends on the rate Xo itself through \i and a. 
Rewritten as 



Xo 



a 



V2tt 
1 1 



exp 



(r > 0) 



-erf 



(38) 



P (jL O (r > 0) is the probability that the intensity of a neuron is 
above threshold and therefore contributes to the transmission of 
a small fluctuation in the input. A neuron for which r < 0 acts as 
if it was absent. Hence we can treat the network with rectifying 
neurons completely analogous to the case of linear Hawkes pro- 
cesses, but multiply the synaptic weight / or — gj of each neuron 
with P^.air > 0), i.e., the linearized connectivity matrix is 



w = P^fr > 0)J. 



(39) 



Figure 7 shows the agreement of the covariance functions 
obtained from direct simulation of the network of Hawkes pro- 
cesses and the analytical solution (21) with average firing rate 
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FIGURE 7 | Covariance structure in spiking networks corresponds to 
OUP with output noise. (A) Autocovariance obtained by direct simulation 
of the LIF (black), Hawkes (gray), and OUP (light gray) models for excitatory 
(dots) and inhibitory neurons (crosses). (B) Covariance C£x averaged over 
disjoint pairs of neurons for LIF (black dots), Hawkes (gray dots), and OUP 
with output noise (empty circles). (C) Covariance averaged over disjoint pairs 
of neurons of the same type. (D) Autocovariance of the population averaged 
activity. Averages in (CD) over excitatory neurons as black dots, over 
inhibitory neurons as gray dots. Corresponding theoretical predictions (21) 
are plotted as light gray curves in all panels except (A). Light gray diagonal 
crosses in (A,D) denote theoretical peak positions determined by the firing 
rate f as rAf (where At = 0.1 ms is the time resolution of the histogram). 



tl the neuron emits an action potential and the membrane poten- 
tial is reset to V r , where it is clamped for the refractory time x r . 
The spiking activity of neuron is described by this sequence 
of action potentials, the spike train s,(f) = £^.8(f — tV). The 
dynamics of a single neuron is deterministic, but in network states 
of asynchronous, irregular activity and in the presence of exter- 
nal Poisson inputs to the network, the summed input to each cell 
can well be approximated as white noise (Brunei, 2000) with first 
moment u,, = x m ^jJijTj and second moment a 2 = x m ^jlffi, 
where r; is the stationary firing rate of neuron The station- 
ary firing rate of neuron i is then given by Fourcaud and Brunei 
(2002) 



r 1 = x r + x m yfn (F{yo) - F(y r )) 

/y 



(41) 



. , V e , r - \Li a x s 

withye.r = h - J ~ 

ai 2 V x m 



with Riemann's zeta function The response of the LIF neuron 
to the injection of an additional spike into afferent j determines 
the impulse response W{jh(t) of the system. The time integral 
Wij = Wy/g 00 h(t) dt is the DC-susceptibility, which can formally 
be written as the derivative of the stationary firing rate by the rate 
of the afferent r ; -, which, evaluated by help of (41), yields (Helias 
et al, 2013, Results and App. A) 



Xo determined by (38), setting the effective strength of the noise 
p 2 = Xo, and the linearized coupling as described above. The 
detailed procedure for choosing the parameters in the direct 
simulation is described together with the implementation of 
the Hawkes model in "Implementation of Hawkes neurons in a 
spiking simulator code". 

5. LEAKY INTEGRATE-AND-FIRE NEURONS 

In this section we consider a network of LIF model neurons with 
exponentially decaying postsynaptic currents and show its equiv- 
alence to the network of OUP with output noise, valid in the 
asynchronous irregular regime. A spike sent by neuron j at time 
t arrives at the target neuron i after the synaptic delay d, elicits a 
synaptic current Z; that decays with time constant x s and causes 
a response in the membrane potential V; proportional to the 
synaptic efficacy / y -. With the time constant x m of the membrane 
potential, the coupled set of differential equations governing the 
subthreshold dynamics of a single neuron i is (Fourcaud and 
Brunei, 2002) 

dV; 

r m -r- = -Vi + hlt) 
dt 

dl N 

where the membrane resistance was absorbed into the defini- 
tions of J,j and I;. If Vj reaches the threshold Ve at time point 



= «/,; + P/f (42) 

with a = ^7t(x m r,) 2 — (/(ye) -f(y r )) 

andp = v^(w,) 2 -^ f/Oe) -f(y r ) . 

2af \ a, a, ) 

In the strongly fluctuation-driven regime, the temporal behavior 
of the kernel h is dominated by a single exponential decay, whose 
time constant can be determined empirically. In a homogeneous 
random network the firing rates of all neurons are identical r; = r 
and follow from the numerical solution of the self-consistency 
Equation (41). Approximating the autocovariance function of a 
single spike train by a 8-peak scaled by the rate fh(t), one obtains 
for the covariance function c between pairs of spike trains the 
same convolution Equation (36) as for Hawkes neurons (Helias 
et al, 2013, cf. equation 5). As shown in "Convolution equation 
for linear noisy rate neurons" this convolution equation coincides 
with that of a linear rate model with output noise (37), where the 
diagonal elements of D are chosen to agree to the average spike 
rate p 2 = r. The good agreement of the analytical cross covari- 
ance functions (21) for the OUP with output noise and direct 
simulation results for LIF are shown in Figure 7. 

6. DISCUSSION 

In this work we describe the path to a unified theoretical view on 
pairwise correlations in recurrent networks. We consider binary 
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FIGURE 8 | Different echo terms for spiking and non-spiking neurons. 

Binary non-spiking neurons shown in (A,C) and LIF in (B,D). (A),(B) Echo 
terms by direct influence of the neuron's output on the network in 
dependence of neuron types (in A,B cee.cei, and Cxx are plotted as 
black, gray dots and circles). (C,D), Contributions to the covariance evoked 
by correlated and common input (black dots) measured with help of 
auxiliary model neurons which do not provide feedback to the network. 
Corresponding theoretical predictions (16) are plotted as light gray curves 
throughout. 



neuron models, LIF models, and linear point process models. 
These models containing a non-linearity (spiking threshold in 
spiking models, non-linear sigmoidal gain function in binary 
neurons, strictly positive rates in Hawkes processes) are lin- 
earized, taking into account the distribution of the fluctuating 
input. 

The work presents results for several neuron models: We derive 
analytical expressions for delay-coupled OUP with input and with 
output noise, we extend the analytical treatment for stochas- 
tic binary neurons to the presence of synaptic delays, present a 
method that takes into account network-generated noise to deter- 
mine the effective gain function, extend the theory of Hawkes 
processes to the existence of delays and inhibition, and present 
in Equation (12) a condition for the onset of global oscillations 
caused by delayed feedback, generalized to feedback pathways 
through different eigenvalues of the connectivity. 

Some results qualitatively extend the existing theory (delays, 
inhibition), others improve the accuracy of existing theories 
(linearization including fluctuations). More importantly, our 
approach enables us to demonstrate the equivalence of each of 
these models after linear approximation to a linear model with 
fluctuating continuous variables. The fact that linear perturba- 
tion theory leads to effective linear equations is of course not 
surprising, but the analytical procedure firstly enables a map- 
ping between models that conserves quantitative results and 
secondly allows us to uncover common structures underlying 
the emergence of correlated activity in recurrent networks. For 
the commonly appearing exponentially decaying response ker- 
nel function, these rate models coincide with the OUP (OUP, 
Uhlenbeck and Ornstein, 1930; Risken, 1996). We find that the 
considered models form two groups, which, in linear approxima- 
tion merely differ by a matrix valued factor scaling the noise and 
in the choice of variables interpreted as neural activity. The differ- 
ence between these two groups corresponds to the location of the 
noise: spiking models — LIF models and Hawkes models — belong 
to the class with noise on the output side, added to the activity of 
each neuron. The non-spiking binary neuron model corresponds 
to an OUP where the noise is added on the input side of each 
neuron. The closed solution for the correlation structure of OUP 
holds for both classes. 

We identify different contributions to correlations in recurrent 
networks: the solution for output noise is split into three terms 
corresponding to the 8-peak in the autocovariance, the covari- 
ance caused by shared input, and the direct synaptic influence of 
stochastic fluctuations of one neuron on another-the latter echo 
terms are equal to propagators acting with delays (Helias et al., 
2013). A similar splitting into echo and correlated input terms for 
the case of input noise is shown in Figure 8. For increasing net- 
work size N — >■ oo, keeping the connection probability p fixed, so 
that K = pN, and with rescaled synaptic amplitudes / ~ l/*/N 
(van Vreeswijk and Sompolinsky, 1996; Renart et al., 2010) the 
echo terms vanish fastest. Formally this can be seen from (18): the 
multiplicative factor of the common covariance term 94 does not 
change with N while the other coefficients decrease. So ultimately 
all four entries of the matrix c have the same time dependence 
determined by the common covariance term 94. In particular 
the covariance between excitation and inhibition csi becomes 



symmetric in this limit. This finally provides a quantitative expla- 
nation of the observation made in (Renart et al, 2010) that the 
time-lag between excitation and inhibition vanishes in the limit 
of infinitely large networks. For a different synaptic rescaling / ~ 
N~ 1 while keeping p 2 constant by appropriate additional input to 
each neuron (see Helias et al, 2013 applied to the LIF model), all 
multiplicative factors decrease ~ N^ 1 and so does the amplitude 
of all covariances. Hence the asymmetry of cgx does not vanish in 
this limit. The same results hold for the case of output noise where 
the term with (pi describes the common input part of the covari- 
ance. In this case and for finite network size, cx£ coincides with 
C££ and C£j with c%x for f > 0, having a discontinuous jump at 
the time of the synaptic delay f = d. For time lags smaller than 
the delay all four covariances coincide. This is due to causality, as 
the second neuron cannot feel the influence of a fluctuation that 
happened in the first neuron less than one synaptic delay before. 
The covariance functions for systems corresponding to an OUP 
with input noise contain neither discontinuities nor sharp peaks 
at t = d, but C£x and c%s have maxima and minima near this 
location. This observation can be interpreted as a result of the 
stochastic nature of the binary model where changes in the input 
influence the state of the neuron only with a certain probability. 
So, the entries of c in this case take different values for |f| < d 
but show the tendency to approach each other with increasing 
\t\ J£> d. This tendency increases with network size. Our analyti- 
cal solutions (18) for input noise and (21) for output noise hence 
explain the model-class dependent differences in the shape of 
covariance functions. 

The two above mentioned synaptic scaling procedures are 
commonly termed "strong coupling" (/ ~ 1 and "weak 
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coupling" (/ ~ 1/N), respectively. The results shown in Figure 6 
were obtained for / = and p" = 0.5, so the number of 

synapses required to cause a notable effect on the gain func- 
tion is 1/(P7) = ^/N, which is small compared to the number 
of incoming synapses pN. Hence the network is in the strong 
coupling regime. Also note that for infinite slope of the gain 
function, p — >■ oo, the magnitude of the covariance becomes 
independent of the synaptic amplitude /, in agreement with the 
linear theory presented here. This finding can readily be under- 
stood by the linearization procedure, presented in the current 
work, that takes into account the network- generated fluctua- 
tions of the total input. The amplitude cr of these fluctuations 
scales linearly in / and the effective susceptibility depends on J /a 
in the case p 1 — »■ oo, explaining the invariance (Grytskyy et al., 
2013). In the current manuscript we generalized this procedure 
to finite slopes P and to other models than the binary neuron 
model. 

Our approach enables us to map results obtained for one neu- 
ron model to another, in particular we extend the theory of all 
considered models to capture synaptic conduction delays, and 
devise a simpler way to obtain solutions for systems considered 
earlier (Ginzburg and Sompolinsky, 1994). Our derivation of 
covariances in spiking networks does not rely on the advanced 
Wiener-Hopf method (Hazewinkel, 2002), as earlier derivations 
(Hawkes, 1971; Helias et al., 2013) do, but only employs elemen- 
tary methods. Our results are applicable for general connectivity 
matrices, and for the purpose of comparison with simulations we 
explicitly derive population averaged results. The averages of the 
dynamics of the linear rate model equations are exact for random 
network architectures with fixed out-degree, and approximate for 
fixed in-degree. Still, for non-linear models the linearization for 



fixed in-degree networks are simpler, because the homogeneous 
input statistics results in an identical linear response kernel for 
all cells. Finally we show that the oscillatory properties of net- 
works of integrate-and-fire models (Brunei, 2000; Helias et al., 
2013) are model-invariant features of all of the studied dynamics, 
given inhibition acts with a synaptic delay. We relate the collective 
oscillations to the pole structure of the cross spectrum, which also 
determines the power spectra of population signals such as EEG, 
ECoG, and the LFP. 

The presented results provide a further step to understand 
the shape and to unify the description of correlations in recur- 
rent networks. We hope that our analytical results will be useful 
to constrain the inverse problem of determining the synaptic 
connectivity given the correlation structure of neurophysiologi- 
cal activity measurements. Moreover the explicit expressions for 
covariance functions in the time domain are a necessary pre- 
requisite to understand the evolution of synaptic amplitudes in 
systems with spike-timing dependent plasticity and extend the 
existing methods (Burkitt et al, 2007; Gilson et al, 2009, 2010) to 
networks including inhibitory neurons and synaptic conduction 
delays. 
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APPENDIX 

CALCULATION OF THE POPULATION AVERAGED CROSS COVARIANCE 
IN TIME DOMAIN 

We obtain the population averaged cross spectrum for the 
Ornstein-Uhlenbeck process with input noise by inserting the 
averaged connectivity matrix w = M (14) into (8). The two eigen- 
values of M are 0 and L = Kw(l — yg). Taking these into account, 
we first rewrite the term 



+ /(.<o)(] +A'.r( Y * _ Y f )i-i,;(o»)MD 



+ DM'/(-ft))(l +Kw 



yg i 
-yg -i 



H d (-co)) + D 



: /(w)MDM T /(-co) +/(w)MD + DM T /(-w) + D, 



(H^y 1 - Mr 1 



Kw Hd(u>) — Kw 



= ((H,j(<o)- J - 0)(H d { M r l - L)) 1 H^oo)" 1 ! + Kw 



Yg ~Yg 
1 -1 



f(w) 1 + Kw 



Yg 



H d (t») 



where we introduced/(co) = (H^co) -1 — L) _1 . The correspond- 
ing transposed and conjugate complex term follows analogously. 
Hence we obtain the expression for the cross spectrum (17). The 
residue off (co) at co = z k (L) is 



Ml — co 



Res(/, a) = z k (L)) = lim — - — - 

(Ol-xo j 1 (CjOi) 



l'Hopital ,. 1 I d(e lu,d (l + iwT)) 

lim 



a>i-*a> (f l )'((x>i) \ dw 

= (ide md (l + icot) + ite" 0 ^" 1 

= (idL + ixj^Y 1 , 

where in the last step we used the condition for a pole 
Hd(zt) 1 = e' Zkd (l + iz k x) = L (see "Spectrum of the dynam- 
ics"). The residue of Hd(w) at z(0) = ^ is — ^e d l x . Using the 
residue theorem, we need to sum over all poles within the inte- 
gration contour {z k (L)\k e N} U l - to get the expression for c(t) = 

i Res(C(z),z)e' zf forf ^ 



/^C(co)rtco = ; £ 



z<E{z k (L)\k€N)U{ 



1 f+OO 

0. Sorting (17) to obtain four matrix prefactors and remain- 
ders with different frequency dependence, OiCco) =/(co)/(— co), 
<J> 2 (co) =/(co)/(— a))Hrf(co), $3(00) = <J>i(— «)> and $4(00) = 
/(«)/(— u>)Hd{ui)Hd(— eo), we get (18). C(co) for output noise 
(20) is obtained by multiplying the expression for C(a>) for input 
noise with H7 (a^HT^— co) = (1 + m 2 t, 2 ). In order to perform 
the back Fourier transformation one first needs to rewrite the 
cross spectrum in order to isolate the frequency independent 
term and the two terms that vanish for either t < d or f > d, as 
described in "Fourier back transformation," 

C(co) = /(w)(l + Kw ~ Y M Hd(w))MDM T /(-w) 

(l+Kw( Yg 1 )H rf (-a>)) 

V-yg -1/ 



where in the last step we used 



M = 0, because M is 



yg -yg 

1 -1 

symmetric, obtaining (21). For each of the first three terms in the 
last expression the right integration contour needs to be chosen 
as described in "Fourier back transformation" on the example of 
the general expression (16). 

IMPLEMENTATION OF NOISY RATE MODELS 

The dynamics is propagated in time steps of duration At (note 
that in other works we use h as a symbol for the computation 
step size, which here is used as the symbol for the kernel). The 
product of the connectivity matrix with the vector of output vari- 
ables at the end of the previous step i — 1 is the vector I(f,) of 
inputs at the current step i. The intrinsic time scale of the system 
is determined by the time constant x. For sufficiently small time 
steps At <SC x these inputs can be assumed to be time independent 
within one step. So we can use (3) or (4) and analytically convolve 
the kernel function h assuming the input to be constant over the 
time interval Af. This corresponds to the method of exponential 
integration (Rotter and Diesmann, 1999, see App. C.6) requir- 
ing only local knowledge of the connectivity matrix w. Note that 
this procedure becomes exact for Af — >■ 0 and for finite A f is an 
approximation. The propagation of the initial value r,-(£j_i) until 
the end of the time interval takes the form r;(tj_i) e~ At / x because 
h(ti) = fe(fj_i) e~ Af//T , so weobtain the expression r;(t,) at the end 
of the step as 



rjifi) 



1) + (1 - e 



-Af/T 



)Ij(ti), 



(43) 



where Ij denotes the input to the neuron j. For output noise the 
output variable of neuron j is y; = rj + x;, with the locally gener- 
ated additive noise Xj and hence the input is L(fi) = (wy(t;))j. In 
the case of input noise the output variable is r; and the additional 
noise is added to the input variable, L(tf) = (wr(f,)) ; + Xj(t{). In 
both cases Xj is implemented as a binary noise: in each time step, Xt 
is independently and randomly chosen to be 1 or —1 with prob- 
ability 0.5 multiplied with p/VA? to satisfy (2) for discretized 
time. Here the 8-function is replaced by a "rectangle" function 
that is constant on the interval of length At, vanishes elsewhere, 
and has unit integral. The factor At -1 in the expression for 
x 2 ensures the integral to be unity. So far, the implementation 
assumes the synaptic delay to be zero. To implement a non-zero 
synaptic delay d, each object representing a neuron contains an 
array b of length Id = d/At acting as a ring buffer. The input Ij(t,) 
used to calculate the output rate at step i according to (43) is then 
taken from position i mod 1^ of this array and after that replaced 
by the input presently received from the network, so that the new 
input will be used only after one delay has passed. This sequence 
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of buffer handling can be represented as 
Ij(ti) <r- b[i mod k] 

(Wl)j + Xj 



b[imoAld\ •«- 



(wy)j 



for input noise 
for output noise 



The model is implemented in Python version 2.7 (Python 
Software Foundation, 2008) using numpy 1.6.1 (Ascher et al., 
2001) and scipy 0.9.0 (Jones et al, 2001). 

IMPLEMENTATION OF BINARY NEURONS IN A SPIKING 
SIMULATOR CODE 

The binary neuron model is implemented in the NEST sim- 
ulator, version 2.2.1 (Gewaltig and Diesmann, 2007), which 
allows distributed simulation on parallel machines and han- 
dles synaptic delays in the established framework for spiking 
neurons (Morrison et al., 2005). The name of the model is 
"ginzburg_neuron". In NEST information is transmitted in 
form of point events, which in case of binary neurons are sent 
if the state of the neuron changes: one spike is sent for a down- 
transition and two spikes at the same time for an up-transition, 
so the multiplicity reflects the type of event. The logic to decode 
the original transitions is implemented in the function handle 
shown in Alg. 2. If a single spike is received, the synaptic weight 
w is subtracted from the input buffer at the position determined 
by the time point of the transition and the synaptic delay. In 
distributed simulations a single spike with multiplicity 2 sent to 
another machine is handled on the receiving side as two sepa- 
rate events with multiplicity 1 each. In order to decode this case 
on the receiving machine we memorize the time (fi ast ) and ori- 
gin (global id gidj ast of the sending neuron) of the last arrived 
spike. If both coincide to the spike under consideration, the send- 
ing neuron has performed an up transition 0 — >■ 1. We hence 
add twice the synaptic weight 2w to the input buffer of the tar- 
get neuron, one that reflects the real change of the system state 
and another that compensates the subtraction of w after recep- 
tion of the first spike of a pair. The algorithm relies on the fact 
that within NEST two spikes that are generated by one neuron at 
the same time point are delivered sequentially to the target neu- 
rons. This is assured, because neurons are updated one by one: 
The update propagates each neuron by a time step equal to the 
minimal delay c? m ; n in the network. All spikes generated within 
one update step are written sequentially into the communication 
buffers, and finally the buffers are shipped to the other processors 
(Morrison et al., 2005). Hence a pair of spikes generated by one 
neuron within a single update step will be delivered consecutively 
and will not be interspersed by spikes from other neurons with 
the same time stamp. 

The model exhibits stochastic transitions (at random points 
in time) between two states. The transitions are governed by 
probabilities $(h). Using asynchronous update (Rumelhart et al., 
1986), in each infinitesimal interval [f, t + of) each neuron in 
the network has the probability ^8f to be chosen for update 
(Hopfield, 1982). A mathematically equivalent formulation draws 
the time points of update independently for all neurons. For a 



particular neuron, the sequence of update points has exponen- 
tially distributed intervals with mean duration x, i.e., it forms a 
Poisson process with rate t~ 1 . We employ the latter formulation 
to incorporate binary neuron models in the globally time-driven 
spiking simulator NEST (Gewaltig and Diesmann, 2007) and con- 
strain the points of transition to a discrete time grid Af = 0. 1 ms 
covering the interval <i m ; n > Af. This neuron state update is 
implemented by the algorithm shown in Alg. 1. Note that the 
field h is updated in steps of Af while the activity state is updated 
only when the current time exceeds the next potential transition 
point. As the last step of the activity update we draw an expo- 
nentially distributed time interval to determine the new potential 
transition time. The potential transition time is represented with 
a higher resolution (on the order of microseconds) than Af to 
avoid a systematic bias of the mean inter-update-interval. This 
update scheme is identical to the one used in (Hopfield, 1982). 
Note that the implementation is different from the classical asyn- 
chronous update scheme (van Vreeswijk and Sompolinsky, 1998), 
where in each discrete time step Af exactly one neuron is picked 
at random. The mean inter-update-interval (time constant x in 
Alg. 1) in the latter scheme is determined by x = AtN, with N 
the number of neurons in the network. For small time steps both 
schemes converge so that update times follow a Poisson process. 

At each update time point the neuron state becomes 1 with 
the probability given by the function 4> applied to the input at 
that time according to (25) and 0 with probability 1 — (f>. The 
input is a function of the whole system state and is constant 
between spikes which indicate state changes. Each neuron there- 
fore maintains a state variable h at each point in time holding 
the summed input and being updated by adding and subtract- 
ing the input read from the ring buffer b at the point readpos(t) 
corresponding to the current time (see Morrison et al, 2005, 
for the implementation of the ring buffer, i.p. Fig 6). The ring 
buffer enables us to implement synaptic delays. For technical 



Algorithm 1 | Update function of a binary neuron embedded in the 
spiking network simulator NEST. 



2 
3 
4 
5 
6 
7 
8 
9 
ID 
11 
12 
13 
14 
15 
16 
17 
18 
111 
20 
21 
22 
23 
24 
25 



y <— 0 // initially neuron is inactive 

r ncxt < r logfrandf)) // next time point of update 

for each time step t: 

h<-h+ f»[readpos(t)] 

if (>«„«,: 

// up-state with probability given by 

// gain function <j> depending on input h(t) 

if (p(h) >rand(): 

else : 

i... «- o 

if !/„,„/!/: 

// down transition : send single spike 
// up transition : send two spikes 
send (ynew +1 spikes ) 



// add an exponentially distributed time interval 
Jiioxt <~ Jnext - rlog(rand()) 



The function readpos(t) returns a position in the ring buffer b corresponding to 
the current time point. 
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Algorithm 2 | Input spike handler of a binary neuron embedded in the 
spiking network simulator NEST. 



1 handlc(i spik „,rf,gid,m) 



2 

3 if m = l: 
4 

5 // mill t ipli city = 1, either a single 1 — > 0 event 

6 // or the first or second of a pair of 0 — > 1 events 
7 

8 i f gid = gid,„, lspikc and f splke = fi„„p Jke : 
9 

10 // received twice the same event , so transition 0 — > 1 

11 // add 2w to compensate for subtraction after reception 

12 // of first event 

13 6[pos(f splke , d. t)\ <- 6[pos(i, piko , d. t)\ + 2w 
14 

15 else: 

16 // count this event negatively , 

17 // assuming it comes as single event 

18 // transition 1 -» 0 

19 6[pos(f, pil «. d, /.)] «- IHI^,, d, ()] - III 

20 else : 
2f 

22 // multiplicity != 1 

23 

24 if in = 2 

25 

26 // count thfs event positively , transition 0 — > 1 

27 b]pas(t, pit ., d., t)\ t- b[pos(t, pik „ d, ()] + w 

28 gMluKpita <- gid 

29 *lastspike *~ ' ; spi k e 



The simulation kernel calls the handle function for each spike event to be deliv- 
ered to the neuron. A spike event is characterized by the time point of occurrence 
fspike. the synaptic delay d after which the event should reach the target, the 
global id gid identifying the sending neuron, and the multiplicity m > 1, indicat- 
ing the reception of multiple spike events. The function pos(f spike , d, t) returns 
the position in the ring buffer b to which the spike is added so that it will be read 
at time t + d by the update function of the neuron, see Alg. 1. 

reasons this implementation requires a minimal delay of a single 
simulation time step (Morrison and Diesmann, 2008). The gain 
function (j) applied to the input h has the form 

$(h) = cih + c 2 - (1 + tanh(c 3 (h - 6))), (44) 

where throughout this manuscript we used cj = 0, ci = 1, and 
c 3 = p, as defined in "Parameters of simulations". 

IMPLEMENTATION OF HAWKES NEURONS IN A SPIKING 
SIMULATOR CODE 

Hawkes neurons (Hawkes, 1971) were introduced in the NEST 
simulator in version 2.2.0 (Gewaltig and Diesmann, 2007). The 
name of the model is "pp_psc_delta". In the following we 
describe the implemented neuron model in general and men- 
tion the particular choices of parameter and correspondences to 
the theory presented in "Hawkes processes". The dynamics of the 
quasi-membrane potential u is integrated exactly within a time 
step Af of the simulation (Rotter and Diesmann, 1999), express- 
ing the voltage «(f,-) at the end of time step i by the membrane 
potential at the end of the previous time step u(f,-i) as 

u(fd = e- At/x Hfe_!) + (1 - e- At/x ) R m I e + b{ti), (45) 

where I e is a time-step wise constant input current (equal to 0 
in all simulations presented in this article) and R m = x m /C m is 
the membrane resistance. The buffer fc(t;) contains the summed 



contributions of incoming spikes, multiplied by their respective 
synaptic weight, which have arrived at the neuron within the 
interval (t,_i, fj]. b is implemented as a ring-buffer in order to 
handle the synaptic delay, logically similar as in "Implementation 
of noisy rate models," described in detail in Morrison et al. (2005). 
The instantaneous spike emission rate is "k = \c\u + C2e C3 "]+, 
where we use C3 = 0 in all simulations presented here. The quan- 
tities in the theory "Hawkes processes," in particular in (35), are 
related to the parameters of the simulated model in the following 
way. The quantity r relates to the membrane potential u as r = 
C\ u + C2 and the background rate v agrees to C2 = v. Hence the 
synaptic weight corresponds to the synaptic weight in the sim- 
ulation multiplied by c\ . For the correspondence of the Hawkes 
model to the OUP with output noise of variance p 2 we use (38) 
to adjust the background rate v in order to obtain the desired rate 
Xrj = p 2 and we choose the synaptic weight / of the Hawkes model 
so that the linear coupling strength w of the OUP agrees to the 
effective linear weight given by (39). These two constraints can 
be fulfilled simultaneously by solving (38) and (39) by numeri- 
cal iteration. The spike emission of the model is realized either 
with or without dead time. In this article we only used the latter. 
In the presence of a dead time, which is constrained to be larger 
than the simulation time step, at most one spike can be generated 
within a time step. A spike is hence emitted with the probability 
p>i = 1 — e XAt , where e XAf is the probability of the comple- 
mentary event (emitting 0 spikes), implemented by comparing 
a uniformly distributed random number to p> 1 . The refractory 
period is handled as described in Morrison et al. (2005). Without 
refractoriness, the number of emitted spikes is drawn from a 
Poisson distribution with parameter XAf, implemented in the 
GNU Scientific Library (Galassi et al, 2006). Reproducibility of 
the random sequences for different numbers of processes and 
threads is ensured by the concept of random number generators 
assigned to virtual processes, as described in (Plesser et al., 2007). 

PARAMETERS OF SIMULATIONS 

For all simulations we used y = 0.25 corresponding to the bio- 
logically realistic fraction of inhibitory neurons, a connectivity 
probability p = 0.1, and a simulation time step of Af = 0.1 ms. 
For binary neurons we measured the covariance functions with a 
resolution of 1 ms, for all other models the resolution is 0.1 ms. 
Simulation time is 10, 000 ms for linear rate and for LIF neurons, 
50, 000 ms for Hawkes, and 100, 000 ms for binary neurons. The 
covariance is obtained for a time window of ±100 ms. 

The parameters for simulations of the LIF model presented 
in Figure 7 and Figure 8 are / = 0.1 mV, x = 20 ms, x s = 2 ms, 
x r = 2 ms, V e = 15 mV, V r = 0, g = 6, d = 3 ms, N = 8000. 
The number of neurons in the corresponding networks of other 
models is the same. Cross covariances are measured between the 
summed spike trains of two disjoint populations of iV rec = 1000 
neurons each. The single neuron autocovariances a a are aver- 
aged over a subpopulation of 100 neurons. The autocovariances 
of the population averaged activity ^-a a + C aa for population 
a G {£ ,2} (shown in Figure 7) are constructed from the esti- 
mated single neuron population averaged autocovariances a a and 
cross covariances C aa . This enables us to estimate a a and C aa 
from the activity of a small subpopulation and still assigns the 
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correct relative weights to both contributions. The correspond- 
ing effective parameters describing the system dynamics are u, = 
15 mV, cr = 10 mV, r = 23.6 Hz (see (40) and the following text 
for details). 

The parameters of the Hawkes model and of the noisy rate 
model with output noise yielding quantitatively agreeing covari- 
ance functions are: 

• For simulations of the noisy rate model with output noise 
presented in Figure 7 and Figure 2 the parameters are w = 
0.0043, g w 5.93, t = 4.07 ms, p 2 = 23.6 Hz, d = 3 ms (see 
(3), (4)). In Figure 2 also results for d = 1 ms and for input 
noise are shown. Signals are measured from N rec = 500 neu- 
rons in each population to obtain cgj, cxs and from the whole 
population to determine egg and cxx- The cross covariances 
Css and Crx are estimated from two disjoint subpopula- 
tions each comprising half of the neurons of the respective 
population. 

• For the network of Hawkes neurons presented in Figure 7 we 
used X 0 w 22.54 Hz (see (38)),/ = 0.0055 mV, d = 3 ms, and 
the same g and x as for the noisy rate model. We measured 
the cross covariances in the same way as for the LIF model, 
but using the spike trains from sub-populations of N rec = 
2000 neurons. The autocovariances of the population averaged 
activity were estimated from the whole populations. 



The network of binary neurons shown in Figure 8 uses 
6 = -3.89mV, P = 0.5mV _1 , / = 0.02mV, d = 3 ms (see 
(25), (44)), and the same g and x as the noisy rate 
model. Covariances are measured using the signals from all 
neurons. 

The simulation results for the network of binary neu- 
rons presented in Figure 6 uses 9 = —2.5 mV, x = 10 ms, P = 
0.5 mV -1 , g = 6, J w 0.0447 mV, N = 2000 and the smallest 
possible value of synaptic delay is d = 0. 1 ms equal to time resolu- 
tion (the same set of parameters only with modified p" = 1 mV" 1 
was used to create Figure 5). The cross covariances Css and Cxi 
are estimated from two disjoint subpopulations each comprising 
half of the neurons of the respective population, csi is measured 
between two such subpopulations. For css and cxx we used the 
full populations. 

The parameters required for a quantitative agreement with 
the rate model with input noise are w « 0.011, p w 2.23 ^ms. 
We used the same parameters in Figure 3, where additionally 
results for w = 0.018 are shown. The population sizes are the 
same as for the binary network. The covariances are estimated 
in the same way as for the rate model with output noise. Note 
that the definition of noisy rate models has no limitation for 
units of p 2 . These can be arbitrary and are chosen differently 
as required by the correspondence with either spiking or binary 
neurons. 
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